/*
 *  pointrsurface.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 17/03/11
 *	Modified by Roberto Ortega on 21/03/11
 *
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "pointrsurface.h"
 
PointRSurface::PointRSurface()
:	Constraint()
{}

// ***
PointRSurface::PointRSurface(vector<string>& vp,
							 int di,
							 vector<double>& lc,
							 Joint* pJoint)
:Constraint(vp,di)
{	
	isContact = true;

	parentJoint = pJoint;
	
	ex = 1.5;
	cs = pJoint->cstiffness;
	cd = pJoint->cdamping;
	cf = pJoint->cfriction;
	
	// Type: STL
	vLocalPoint = lc;
	
	//-- Bisection Method --//
	tole = 1E-10;
	max = 100;

	//-- Bodies --//
	Body1 = pJoint->PBody1;
	Body2 = pJoint->PBody2;
	
	//-- Surface of BodyB --//
	surface = pJoint->PBody2->surface;
	
	//-- Contact --//
	V_old = 0.;
	S_old = 0.;
	
	V = V_old;
	S = S_old;
}

// ***
void PointRSurface::calcFqC(const int& nt)
{
	int i;
	int j;
	int k;

	Fq.setZero(numDofs);
	
	if(integratorType == 2)
	{
		dqm = 0.5 * (dq + dqC);
	}
	else
	{
		dqm = dq;
	}	
	
	Fs.setZero();
	Fd.setZero();
	Ff.setZero();
	
	if(integratorType == 2)
	{
		calV(q);
		
		DV = V - V_old;
		DS = S - S_old;

		//-- Elastic Force --//
		if (DV != 0.0 and DS != 0.0 and cs != 0.0)
		{
			computeForceStiffness();
		}	
		
		//-- Damping Force --//
		if (DV != 0.0 and DS != 0.0 and cs != 0.0 and cd != 0.0)
		{
			computeForceDamping();
		}
		
		//-- Friction Force --//
		if (DV != 0.0 and DS != 0.0 and cs != 0.0 and cf != 0.0)
		{
			computeForceFriction();
		}
	}
	
	if(integratorType != 2)
	{
		calV(q);

		//-- Elastic Force --//
		if (V != 0.0 and cs != 0.0)
		{
			computeForceStiffness();
		}
		
		//-- Damping Force --//
		if (V != 0.0 and cs != 0.0 and cd != 0.0)
		{
			computeForceDamping();
		}
		
		//-- Friction Force --//
		if (V != 0.0 and cs != 0.0 and cf != 0.0)
		{
			computeForceFriction();
		}		
	}
	
	//-- Asemmble all forces --//
	Fq = Fs + Fd + Ff;
}

// ***
void PointRSurface::calcKqC(int nt)
{
	// Using Numerical tangent.
	calcLocalKqNumeric();
	
	// Analytical tangent, TO BE DONE
	/*
	 int i;
	 int j;
	 
	 for (i = 0; i < numDofs; ++i)
	 for (j = 0; j < numDofs; ++j)
	 Kq[i][j] = 0.0;
	 */	
}

// ***

void PointRSurface::calcKdqC(double & partial_qdot)
{
	//cout << "partial_qdot: " << partial_qdot;
}

// ***
void PointRSurface::initializeConsResC(const bool& st, const double& grav, const int& gravDir, const bool augl)
{
	int i;
	
	// numDofs = number of points * dim
	
    Fq.setZero(numDofs);
	qBeta.resize(numDofs);
	
	dqm.resize(numDofs);
	dV.resize(numDofs);
	
	Fs.resize(numDofs);
	Fd.resize(numDofs);
	Ff.resize(numDofs);
	
	// Position vector for the contact point referred to system in A.
	// Local (3D: 3x1) // Type: STL
	aL.setZero();
	for (i = 0; i < 3; ++i)
		aL(i) = vLocalPoint[i];
	
	//Global Position (3D: 3x1)
	aG.setZero();
	bG.setZero();

	// Global Velocity (3D: 3x1)
	daG.setZero();
	dbG.setZero();
	
	// Local Position vector for the contact point referred to system in B.
	// Local (3D: 3x1)
	bL.setZero();
	
	// Base attached to Body1 (qA)
	if (!Body1->type.compare("rigidbody3d4p"))
	{
		// Body1: RigidBody (3D: 4 * 3 = 12)
		qA.resize(4 * dim);
		dqA.resize(4 * dim);
		CA = getC(aL, 4);
	}
	else
	{
		// Body1: FlexiBody (3D: 1 * 3 = 3)
		qA.resize(1 * dim);
		dqA.resize(1 * dim);
		CA = getC(aL, 1);
	}

	J.resize(dim, numDofs);
	D.resize(dim, numDofs);

	// Base attached to Body2 (qB)
	if (!Body2->type.compare("ground") or !Body2->type.compare("rigidbody3d4p"))
	{
		// Body2: Ground or RigidBody (3D: 4 * 3 = 12)
		qB.resize(4 * dim);
		dqB.resize(4 * dim);
		
		C0B = getC(Eigen::Vector3d::Zero(), 4);
		C1B = getC(Eigen::Vector3d::UnitX(), 4);
		C2B = getC(Eigen::Vector3d::UnitY(), 4);
		C3B = getC(Eigen::Vector3d::UnitZ(), 4);
	}
}

// ***

void PointRSurface::initializeConsTangC(const bool& st)
{	
	Kq.setZero(numDofs,numDofs);
}

// ***
void PointRSurface::calcEnerC(double& Ekin, double& Epot, double& Econ, double& Eamb, double& grav, int& gravDir, vector<double>& q_in, vector<double>& dq_in)
{
	Econ += V;
}

//  ***
void PointRSurface::updateCons(void)
{
	//if(Fs.norm() != 0) cout << "||Fs||: " << Fs.norm() << endl;
	//if(Fd.norm() != 0) cout << "||Fd||: " << Fd.norm() << endl;
	//if(Ff.norm() != 0) cout << "||Ff||: " << Ff.norm() << endl << endl;

	if(integratorType == 2 and cs != 0.0)
	{
		V_old = V;
		S_old = S;
	}
	if(integratorType == 2 and cd != 0.0)
	{
		V_old = V;
		S_old = S;
	}
	if(integratorType == 2 and cf != 0.0)
	{
		V_old = V;
		S_old = S;
	}	
}

// ***
void PointRSurface::calcPhi(const VectorXd & Q)
{}

// ***
void PointRSurface::calcDPhi(const VectorXd & Q)
{}

// ***
void PointRSurface::calcDDPhi(const VectorXd & Q)
{}

// ***
void PointRSurface::setFqSymmetries()
{}

// ***
void PointRSurface::setKqSymmetries()
{}


//* Begin Local functions *****************************************************!

// ***
void PointRSurface::makeM()
{}

// ***
void PointRSurface::computeForceStiffness(void)
{
	if (integratorType == 2)
		caldVBeta();

	if (integratorType != 2)
		caldV(q);
	
	Fs = cs * G * dG * dS * J;
}

// ***
void PointRSurface::computeForceDamping(void)
{
	D = (dS * J).transpose() * (dS * J);
	
	//Fd = cd * dG * dG * D * dqm;
	//Fd = cd * dG * dG * D * dqm / (dS * J).norm();
	
	if ((D * dqm).norm() != 0.0)
	{
		//Fd = cd * dG * dG * D * dqm / (D * dqm).norm();
		Fd = cd * dG * dG * D * dqm;
	}
	else
	{
		Fd.setZero();
	}
}

// ***
void PointRSurface::computeForceFriction(void)
{
	int i;
	
	//-- Base attached to Body1 (qA) --//
	if (!Body1->type.compare("rigidbody3d4p"))
	{
		// Body1: RigidBody
		for (i = 0; i < dqA.size(); ++i)
			dqA(i) = dq[i];			
		
		// Get the Global Velocity of Point "a" (Body1 - qA)
		getGVel(dqA, aL, daG);
	}
	else
	{
		// Body1: FlexiBody
		for (i = 0; i < dqA.size(); ++i)
			daG(i) = dq[i];
	}
	
	//-- Base attached to Body2 (qB) --//
	if (!Body2->type.compare("ground"))
	{
		// Body2: Ground
		for (i = 0; i < qB.size(); ++i)
			dqB(i) = 0.0;
		
		// Get the Global Velocity of Point "b" (Body2 - qB)
		getGVel(dqB, bL, dbG);
		
	}
	else
	{
		// Body2: RigidBody
		for (i = 0; i < dqB.size(); ++i)
			dqB(i) = dq[dqA.size() + i];
		
		// Get the Global Velocity of Point "b" (Body2 - qB)
		getGVel(dqB, bL, dbG);
	}
	
	double gamma;
	
	Eigen::Vector3d vt;
	Eigen::Vector3d tt;
	Eigen::Vector3d nn;
	Eigen::Matrix3d ii;
	
	nn = dS / dS.norm();
	ii = Eigen::MatrixXd::Identity(3,3);
	
	vt = (ii - nn * nn.transpose()) * (daG - dbG);
	//vt = (ii - nn * nn.transpose()) * (daG);
	
	if( vt.norm() != 0.0)
	{
		tt = vt / vt.norm();
		//gamma = tanh(vt.norm() / 0.5);
		gamma = vt.norm() / sqrt(vt.norm() * vt.norm() + 0.025);
		//gamma = vt.norm() / sqrt(vt.norm() + 0.1);
	}
	else
	{
		tt.setZero();
		gamma = 0.0;
	}
	
	//Ff = cf * (Fs + Fd).norm() * gamma * tt.transpose() * J;
	//Ff = cf * (Fs + Fd).norm() * gamma * tt.transpose() * J / (tt.transpose() * J).norm();
	
	if ((tt.transpose() * J).norm() != 0.0)
	{
		Ff = cf * (Fs + Fd).norm() * gamma * tt.transpose() * J / (tt.transpose() * J).norm();
		//Ff = cf * (Fs + Fd).norm() * gamma * tt.transpose() * J;
	}
	else
	{
		Ff.setZero();
	}	
}

// ***
void PointRSurface::setCoordinates(VectorXd& q_in)
{
	// q_in (Type: Eigen)
	
	int i;
	
	//-- Base attached to Body1 (qA) --//
	if (!Body1->type.compare("rigidbody3d4p"))
	{
		// Body1: RigidBody
		for (i = 0; i < qA.size(); ++i)
			qA(i) = q_in(i);			
		
		// Get the Global Position of Point "a" (Body1 - qA)
		getGPos(qA, aL, aG);
	}
	else
	{
		// Body1: FlexiBody
		for (i = 0; i < qA.size(); ++i)
			aG(i) = q_in(i);
	}
	//-- Base attached to Body2 (qB) --//
	if (!Body2->type.compare("ground"))
	{
		// Body2: Ground
		for (i = 0; i < qB.size(); ++i)
			qB(i) = Body2->localSystem[i];
	}
	else
	{
		// Body2: RigidBody
		for (i = 0; i < qB.size(); ++i)
			qB(i) = q_in(qA.size() + i);
	}	
	// Get the Local Position of Point "b" (Body2 - qB)
	// Using Global Position of Point "a"
	
	getLPos(qB, bL, aG);
}

// ***
void PointRSurface::calV(VectorXd& q_in)
{

	// q_in (Type: STL)
	setCoordinates(q_in);
	
	// Evaluate Surface in Local Coordinates evalS(xL, yL, zL)
	S = surface->evalS(bL.x(), bL.y(), bL.z());

	// Evaluate Gap
	G = getGap(S);

	// Evaluate Contact Potential V
	V = 0.5 * cs * G * G;
}

// ***
void PointRSurface::caldV(VectorXd & q_in)
{
	// q_in (Type: Eigen)
	
	int i;
	int j;
	
	// q_in (Type: STL)
	setCoordinates(q_in);
	
	// Evaluate Surface and Gradient of Surface in Local Coordinates
	S = surface->evalS(bL.x(), bL.y(), bL.z());
	dS = surface->evaldS(bL.x(), bL.y(), bL.z());
	
	// Evaluate Gap and dG
	G = getGap(S);
	dG = getdGap(S);
	
	RB = getR(qB);
	calMatJ();
	
	dV = cs * G * dG * dS * J;
}

// ***
void PointRSurface::calddV(VectorXd & q_in)
{
	/*
	int i;
	int j;
	int k;
	int l;
	
	setCoordinates(q_in);
	
	// Evaluate Surface and Gradient of Surface in Local Coordinates
	S = surface->evalS(vPosBL[0], vPosBL[1], vPosBL[2]);
	dS = surface->evaldS(vPosBL[0], vPosBL[1], vPosBL[2]);
	ddS = surface->evalddS(vPosBL[0], vPosBL[1], vPosBL[2]);
	
	// Evaluate Gap and dG
	G = getGap(S);
	dG = getdGap(S);
	ddG = getddGap(S);
	
	mRB = getMatR(qB);
	calMatJ();
	
	for (i = 0; i < numDofs; ++i)
		for (j = 0; j < numDofs; ++j)
		{
			ddV[i][j] = 0.0;

			for (k = 0; k < 3; ++k)
			{
				ddV[i][j] += dS[k] * mJ[k][i] * dS[k] * mJ[k][j] ;
				ddV[i][j] *= cs * (dG * dG + G * ddG);
			}
			
			for (k = 0; k < 3; ++k)
				for (l = 0; l < 3; ++l)
				{
					ddV[i][j] += cs * dG * mJ[k][i] * ddS[k][l] * mJ[l][j] ;
				}
		}
	 */
}

void PointRSurface::calMatJ(void)
{
	int j;
	int c;
	
	bG = aG - qB.head(3);	
	
	e1B = RB.col(0);
	e2B = RB.col(1);
	e3B = RB.col(2);
	
	j = 0;
	c = qA.size();
	
	//(1x3)([12,3])
	J.block(0, j, 1, c) = e1B.transpose() * CA;
	J.block(1, j, 1, c) = e2B.transpose() * CA;
	J.block(2, j, 1, c) = e3B.transpose() * CA;
	
	if (!Body2->type.compare("rigidbody3d4p"))
	{
		j = qA.size();
		c = qB.size();
		
		//(1x3)([12,0])
		J.block(0, j, 1, c) = bG.transpose() * C1B - (bG + e1B).transpose() * C0B;
		J.block(1, j, 1, c) = bG.transpose() * C2B - (bG + e2B).transpose() * C0B;
		J.block(2, j, 1, c) = bG.transpose() * C3B - (bG + e3B).transpose() * C0B;
	}
}

// ***
double PointRSurface::getGap(double S)
{
	double Gap;

	if(S < 0.0)
	{
		Gap = pow(-S, ex);
	}
	else
	{
		Gap = 0.0;
	}
	
	return Gap;
}

// ***
double PointRSurface::getdGap(double S)
{
	double dGap;
	
	if(S < 0.0)
	{
		dGap = - ex * pow(-S, ex - 1.0);
	}
	else
	{
		dGap = 0.0;
	}
	
	return dGap;
}

// ***
double PointRSurface::getddGap(double S)
{
	double ddGap;
	
	if(S < 0.0)
	{
		ddGap = ex * (ex - 1.0) * pow(-S, ex - 2.0);
	}
	else
	{
		ddGap = 0.0;
	}
	
	return ddGap;
}

// ***
void PointRSurface::caldVBeta(void)
{
	int ite;
	
	double a;
	double b;
	double c_new;
	double c_old;
	
	double error;
	
	a = 0.0;
	b = 1.0;
	c_new = 0.5;
	c_old = 0.5;
	
	ite = 1;
	do
	{
		c_new = (a + b) / 2.0;
		
		if (fBeta(c_new) == 0.0)
		{
			break;
		}
		if(fBeta(a) * fBeta(c_new) < 0.0)
		{
			b = c_new;
		}
		else
		{
			a = c_new;
		}
		
		if(c_new != c_old)
		{
			error = fabs(c_new - c_old)/ fabs(c_new);
		}
		else
		{
			error = fabs(fBeta(c_new));
		}
		
		c_old = c_new;
		
		ite++;
	}
	while(error > tole && ite < max);
	
	beta = c_new;
	fBeta(beta);
	
	if (ite == max && error > tole)
	{
		cout << "======= Bisection-Beta Method Fails =======" << endl;
		cout << fBeta(0.0) << "\t" << fBeta(1.0) << endl;		
		cout << ite << "\t" << beta << "\t" << fBeta(beta) << endl;
		cout << "======= Bisection-Beta Method Fails =======" << endl;
	}
	if (ite <= max && error < tole)
	{
		//		cout << "========== Bisection Method ==========" << endl;
		//		cout << fBeta(0.0) << "\t" << fBeta(1.0) << endl;
		//		cout << ite << "\t" << beta << "\t" << fBeta(beta) << endl;
		//		cout << "========== Bisection Method ==========" << endl;
	}	
}


// ***
void PointRSurface::calNewtonBeta(void)
{
	int ite;
	
	double fa;
	double fb;
	double ep;
	
	double df;
	double dc;
	
	double xc;
	
	ep = 1E-06;
	xc = 0.5;
	
	ite = 1;
	
	do
	{
		fa = fBeta(xc - ep);
		fb = fBeta(xc + ep);
		
		df = (fb - fa) / (2 * ep);
		
		dc = - fBeta(xc) / df;
		
		xc = xc + dc;
		
		ite++;
	}
	while((fabs(dc) > tole) && (ite < max));
	
	beta = xc;
	fBeta(beta);
	
	if (ite == max && fabs(dc) > tole)
	{
		cout << "======= Newton Method-Beta Fails =======" << endl;
		cout << fBeta(0.0) << "\t" << fBeta(1.0) << endl;		
		cout << ite << "\t" << beta << "\t" << fBeta(beta) << endl;
		cout << "======= Newton Method-Beta Fails =======" << endl;
	}
	if (ite <= max && fabs(dc) < tole)
	{
//		cout << "========== Newton-Beta Method ==========" << endl;
//		cout << fBeta(0.0) << "\t" << fBeta(1.0) << endl;		
//		cout << ite << "\t" << beta << "\t" << fBeta(beta) << endl;
//		cout << "========== Newton-Beta Method ==========" << endl;
	}
}


// ***
double PointRSurface::fBeta(double beta)
{	
	int i;
	double R;

	// qBeta, q, qC (Type: Eigen)
	qBeta = qC + beta* (q - qC);
	
	caldV(qBeta);
	
	R = 0.0;
	for (i = 0; i < numDofs; ++i)
		R += dV(i) * (q(i) - qC(i));
	/*
	 double e = 0.0;
	 
	 if(DV >= 0 )
	 //return (R - DV);
	 return (R - (e - 1) * DV);
	 else
	 return (R - e * DV);
	 */
	
	return (R - DV);
}

//*	Temporally functions ******************************************************!

void PointRSurface::calcLocalKqNumeric()
{
	int i;
	int j;
	
	double eps = perturbation;
	
	//eps = 1E-8;
	
	if(Fq.size()!=0)
	{
		if(Kq.size()==0)
		{
			Kq.setZero(numDofs, numDofs);
		}
		for (j = 0; j < numDofs; ++j)
		{
			q(j) += eps;
		    // 1 means no simultaneous tangent computation (if available)
			calcFqC(1);
			
			for (i = 0; i < numDofs; ++i)
				Kq(i,j) = Fq(i);
			
			q(j)-= 2.0 * eps;
			// 1 means no simultaneous tangent computation (if available)
			calcFqC(1);
			
			for (i = 0; i < numDofs; ++i)
				Kq(i,j) -= Fq(i);
			
			for (i = 0; i < numDofs; ++i)
				Kq(i,j) /= (2.0 * eps);
			
			q(j) += eps;
		}
		// 1 means no simultaneous tangent computation (if available)
		calcFqC(1);
	}
}

//**
Eigen::Matrix3d PointRSurface::getR(Eigen::VectorXd& q)
{
	// Assumed that q contains the coordinates of 4 points that define an
	// orthonormal base
	// First point is the origin of base
	
	Eigen::Matrix3d R;
		
	R(0,0) = q(3) - q(0);
	R(0,1) = q(6) - q(0);
	R(0,2) = q(9) - q(0);
	
	R(1,0) = q(4) - q(1);
	R(1,1) = q(7) - q(1);
	R(1,2) = q(10) - q(1);
	
	R(2,0) = q(5) - q(2);
	R(2,1) = q(8) - q(2);
	R(2,2) = q(11) - q(2);
	
	return R;
}

//**
Eigen::MatrixXd PointRSurface::getC(Eigen::Vector3d v, int size)
{
	int i;
	int j;
	int k;
	
	Eigen::MatrixXd C;
	Eigen::MatrixXd I;
	
	C.setZero(dim, size * dim);
	I.setIdentity(dim, dim);
	
	if (size == 1)
	{
		C.block(0,0,dim,dim) = I * (1 - v.x() - v.y() - v.z());
		 
	}
	
	if (size == 4)
	{
		C.block(0,0,dim,dim) = I * (1 - v.x() - v.y() - v.z());
		C.block(0,3,dim,dim) = I * v.x();
		C.block(0,6,dim,dim) = I * v.y();
		C.block(0,9,dim,dim) = I * v.z();

	}
	
	return C;
}

//**
void PointRSurface::getGPos(Eigen::VectorXd& q,
							Eigen::Vector3d& vL,
							Eigen::Vector3d& vG)
{
	// Assumed that "q" contains the coordinates of 4 points
	// that define an orthonormal base
	// First point is the origin of base
	
	// r = r0 + rho
	// r0 = q(0:2)
	// rho = vLx * (r1 - r0) + vLy * (r2 - r0) + vLz * (r3 - r0)
	
	vG.setZero();
	
	vG.x() += q(0);
	vG.x() += vL.x() * (q(3) - q(0));
	vG.x() += vL.y() * (q(6) - q(0));
	vG.x() += vL.z() * (q(9) - q(0));
	
	vG.y() += q(1);
	vG.y() += vL.x() * (q(4) - q(1));
	vG.y() += vL.y() * (q(7) - q(1));
	vG.y() += vL.z() * (q(10) - q(1));
	
	vG.z() += q(2);
	vG.z() += vL.x() * (q(5) - q(2));
	vG.z() += vL.y() * (q(8) - q(2));
	vG.z() += vL.z() * (q(11) - q(2));	
}

//**
void PointRSurface::getGVel(Eigen::VectorXd& dq,
							Eigen::Vector3d& vL,
							Eigen::Vector3d& dvG)
{
	// Assumed that "dq" contains the velocities of 4 points
	// that define an orthonormal base
	// First point is the origin of base
	
	// dr = dr0 + drho
	// dr0 = dq(0:2)
	// rho = vLx * (dr1 - dr0) + vLy * (dr2 - dr0) + vLz * (dr3 - dr0)
	
	dvG.setZero();
	
	dvG.x() += dq(0);
	dvG.x() += vL.x() * (dq(3) - dq(0));
	dvG.x() += vL.y() * (dq(6) - dq(0));
	dvG.x() += vL.z() * (dq(9) - dq(0));
	
	dvG.y() += dq(1);
	dvG.y() += vL.x() * (dq(4) - dq(1));
	dvG.y() += vL.y() * (dq(7) - dq(1));
	dvG.y() += vL.z() * (dq(10) - dq(1));
	
	dvG.z() += dq(2);
	dvG.z() += vL.x() * (dq(5) - dq(2));
	dvG.z() += vL.y() * (dq(8) - dq(2));
	dvG.z() += vL.z() * (dq(11) - dq(2));	
}

//**
void PointRSurface::getLPos(Eigen::VectorXd& Q,
							Eigen::Vector3d& vL,
							Eigen::Vector3d& vG)
{
	// Assumed that q contains the coordinates of 4 points that define an orthonormal base
	// First point is the origin of base
	// alpha = [ (r-r0) * (r1-r0) ] / (r1-r0)^2  (with (r1-r0)^2=1)
	
	// alpha = (r - r0) * (r1 - r0)
	// beta =  (r - r0) * (r2 - r0)
	// gamma = (r - r0) * (r3 - r0)
	
	vL.setZero();
		
	vL.x() += (vG.x() - Q(0)) * (Q(3) - Q(0));
	vL.x() += (vG.y() - Q(1)) * (Q(4) - Q(1));
	vL.x() += (vG.z() - Q(2)) * (Q(5) - Q(2));
	
	vL.y() += (vG.x() - Q(0)) * (Q(6) - Q(0));
	vL.y() += (vG.y() - Q(1)) * (Q(7) - Q(1));
	vL.y() += (vG.z() - Q(2)) * (Q(8) - Q(2));

	vL.z() += (vG.x() - Q(0)) * (Q(9) - Q(0));
	vL.z() += (vG.y() - Q(1)) * (Q(10) - Q(1));
	vL.z() += (vG.z() - Q(2)) * (Q(11) - Q(2));
}

/***End File ******************************************************************/
